{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 吉布斯采样"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**吉布斯采样**可以被看做Metropolis-Hastings算法的⼀个具体的情形。吉布斯采样的每个步骤涉及到将⼀个变量的值替换为以剩余变量的值为条件，从这个概率分布中抽取的那个变量的值。因此我们将$z_i$替换为从概率分布$p(z_i|z_{\\setminus i})$中抽取的值，这个步骤要么按照某种特定的\n",
    "顺序在变量之间进⾏循环，要么每⼀步中按照某个概率分布随机地选择⼀个变量进⾏更新。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "n维吉布斯采样\n",
    "\n",
    "* 初始化$\\{z_i:i=1,...,M\\}$\n",
    "* 对于$\\tau=1,...,T$:\n",
    "     * 采样$z_1^{\\tau+1}\\sim p(z_1|z_2^{\\tau},z_3^{\\tau},...z_M^{\\tau})$。\n",
    "     * 采样$z_2^{\\tau+1}\\sim p(z_2|z_1^{\\tau+1},z_3^{\\tau},...z_M^{\\tau})$。\n",
    "     * $\\vdots$  \n",
    "     * 采样$z_j^{\\tau+1}\\sim p(z_j|z_1^{\\tau+1},...z_{j-1}^{\\tau+1},z_{j+1}^{\\tau},...,z_M^{\\tau})$。\n",
    "     * $\\vdots$   \n",
    "     * 采样$z_M^{\\tau+1}\\sim p(z_M|z_1^{\\tau+1},z_2^{\\tau+1},...,z_{m-1}^{\\tau+1})$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "吉布斯采样步骤看成Metropolis-Hastings算法的⼀个特定的情况。若⼀个Metropolis-Hastings采样的步骤，它涉及到变量$z_k$，同时保持剩余的变量$z_{\\setminus k}$不变，从$z$到$z^*$的转移概率为$q_k(z^*|z)=p(z_k^*|z_{\\setminus k})$。\n",
    "\n",
    "因为在采样过程中向量的各个元素不变，所以$z^*_{\\setminus k}=z_{\\setminus k}$,又由$p(z)=p(z|z_{\\setminus k})p(z_{\\setminus k})$,那么Metropolis-Hastings算法中的接受概率的因⼦为:\n",
    "$$A(z^*,z)=\\frac{p(z^*)q_k(z|z*)}{p(z)q_k(z^*|z)}=\\frac{p(z_k^*|z^*_{\\setminus k})p(z^*_{\\setminus k})p(z_k|z^*_{\\setminus k})}\n",
    "{p(z_k|z_{\\setminus k})p(z_{\\setminus k})p(z^*|z_{\\setminus k})}=1$$\n",
    "\n",
    "所以Metropolis-Hastings步骤总是被接受的"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy.stats as st\n",
    "\n",
    "def finite_mixture_model():\n",
    "    # Generate some data\n",
    "    k = 3\n",
    "    # target distribution:3 2D-Gaussian distribution\n",
    "    mu = np.array([[2,3],[6,4],[4,5]])\n",
    "    cov = np.array([[[1,0],[0,1]],\n",
    "                   [[1,1],[-1,1]],\n",
    "                   [[1,-1],[0,1]]])\n",
    "    #sampling frequence\n",
    "    freq = [200, 100, 50]\n",
    "    colors = ['ro', 'bo', 'go']\n",
    "\n",
    "    X = np.zeros((0,2))\n",
    "\n",
    "    for i in range(k):\n",
    "        vals = np.random.multivariate_normal(mu[i], cov[i], freq[i])\n",
    "        plt.plot(vals[:,0], vals[:,1], colors[i])\n",
    "        X = np.vstack([X, vals])\n",
    "        print(X.shape)\n",
    "    plt.show()\n",
    "\n",
    "    # Constants\n",
    "    alpha = np.ones(k)/k\n",
    "    sigma = 1\n",
    "    rho = 1\n",
    "    n,d = X.shape\n",
    "\n",
    "    # Initialize variables\n",
    "    pi = np.ones(k)/k\n",
    "    z = np.random.multinomial(1, pi, size=n)\n",
    "    mu = np.zeros((k, d))\n",
    "\n",
    "    # Gibbs Sampler\n",
    "    num_iter = 1000\n",
    "    burn_in = 200\n",
    "    samples = []\n",
    "    for count in range(num_iter):\n",
    "        #print (count)\n",
    "\n",
    "        # Update pi\n",
    "        n_k = np.sum(z.T, axis=1, dtype=float)\n",
    "        pi = np.random.dirichlet(alpha + n_k, 1)\n",
    "\n",
    "        # Update mean and covariance of condotional Gaussian distribution\n",
    "        x_hat = np.dot(z.T, X) / n_k[:, None]\n",
    "        mu_hat = x_hat/(1+sigma/(rho*n_k[:, None]))\n",
    "        cov_hat = sigma*rho/(sigma+rho*n_k)\n",
    "        mu = np.array([np.random.multivariate_normal(mu_hat[i], cov_hat[i]*np.eye(d), 1)[0] \\\n",
    "                       for i in range(k)])\n",
    "\n",
    "        # Update z\n",
    "        z_hat = np.zeros((n,k))\n",
    "        for i in range(n):\n",
    "            for c in range(k):\n",
    "                z_hat[i,c] = np.prod(st.norm.pdf(X[i], mu[c], rho*np.ones(d)))\n",
    "        z_hat = (z_hat*pi[:, None])[0]\n",
    "        for i,x in enumerate(z_hat):\n",
    "            z[i] = np.random.multinomial(1, x/sum(x))\n",
    "\n",
    "        samples.append(mu)\n",
    "\n",
    "    # Remove burn in\n",
    "    samples = np.array(samples[burn_in:])\n",
    "\n",
    "    # Find the median of each mu\n",
    "    medians = np.median(samples, axis=0)\n",
    "    plt.plot(medians[0,0], medians[0,1], 'rs')\n",
    "    plt.plot(medians[1,0], medians[1,1], 'gs')\n",
    "    plt.plot(medians[2,0], medians[2,1], 'bs')\n",
    "    plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(200, 2)\n",
      "(300, 2)\n",
      "(350, 2)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/leilichuan/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:20: RuntimeWarning: covariance is not positive-semidefinite.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnW2MHdd53//PLndrXtIh7Su1SC3vXX8oHLim38S2Tl0Y\nSpmgsWrZdT4YNq5cVi6wMdumZFMjqLWtKRfeFm2DSgRaBVg4ChTdWxuu6sS1qzixFRdIgVop5RfR\ntpI0dbi03LiRqZq0vGpIkU8/zA45d+45M+fMnHm9/x8w2N2783Jm5s7/PPM8z3mOqCoIIYR0h6Wm\nG0AIIcQPCjchhHQMCjchhHQMCjchhHQMCjchhHQMCjchhHQMCjchhHQMCjchhHQMCjchhHSMfVXs\n9JZbbtH19fUqdk0IIb3kySef/J6q3uqybiXCvb6+jrNnz1axa0II6SUisuO6Ll0lhBDSMSjchBDS\nMSjchBDSMSjchBDSMSjchBDSMSjchARmem6K9QfWsfSRJaw/sI7puWnTTSI9o5J0QEIWlem5KTY+\ns4Hdq7sAgJ1LO9j4zAYAYHxk3GTTSI+gxU1IQDYf37wh2jG7V3ex+fhmQy0ifYTCTUhALly64PU5\nIUWgcBMSkLVDa16fE1IECjchAdk6toXBymDms8HKAFvHthpqEekjFG5CAjI+Msb2XdsYHRpBIBgd\nGmH7rm0GJklQRFWD7/To0aPKIlOEEOKOiDypqkdd1qXFTQghHYPCTQghHYPCTQghHYPCTQghHYPC\nTQghHYPCTQghHYPCTQghHSNXuEXk1SLy1cRyWURO1dE4Qggh8+SWdVXV3wfwBgAQkWUA3wHwaxW3\nixBCiAVfV8kxAP9LVZ2nkSeEEBIWX+F+D4CPV9EQQgghbjgLt4isAngHgP9o+f+GiJwVkbPPPvts\nqPYRQghJ4WNxvw3Al1X1/5j+qarbqnpUVY/eeuutYVpHSAuYToH1dWBpKfo55RSSpGF8hPu9oJtk\nIeBktzeZToGNDWBnB1CNfm5sULxJszgJt4gcAPBTAD5VbXNI08ST3e5c2oFCb0x2u6jivbkJ7M5O\nIYnd3ejzuqHlT2KchFtVf6iqQ1W9VHWDSLNwsttZLlimirR9XhW0/EkSjpwkM3Rpsts6XDprlqki\nbZ9XRZssf9I8FG4yQ1cmu63LpbO1BQxmp5DEYBB9XidtsfxJO6Bwkxm6MtltXS6d8RjY3gZGI0Ak\n+rm9HX1eFSZfdlssf9IOKNxkhq5MdhvCpePqahmPgfPngevXo59Vi7bJl33nne2w/Ek74GTBpJOs\nP7COnUvzlRdGh0Y4f+p87vaxqyVptQ9WBo13UuvrkVinGY0ikd7cjNwja2vR31V2IqRefCYLpnCT\nTlJWeG3CvyzLuK7XsXZoDVvHtmoX8aWlyNJOIxJZ/KS/cJZ30nvKunRsLpVreq3R/HX6sokLtLjJ\nQnLLv74FF1+4mLueq+slFLGPO5n6NxhUHxAlzUOLm5AMpuemuPynl53WrTt/vYksFtI9cidSIKRv\nbD6+iavXrzqt20T++nhMoSbZ0OImC4erFd3G/HVCAAo3aZgmKhHarOjh/mHr89cJAegqIQ2STumL\nMzkAVCqYW8e2jKmEZ952hkJNOgEtbtIYTVUi7MroUEJs0OImjdFkJcLxkTGFmnQWWtzESC0lUztS\niZCQtkHhJnPUVjK1I5UICWkbFG4yR20lU+lrJqQQ9HGTOer0PdPXXIzplJUCFxla3GQO+p7bDeef\nJK6zvB8WkUdF5PdE5GkR+fGqG0aag77nWZoYJJQF558krhb3GQCfU9UfA/B6AE9X1yTSNPQ936Su\nQO2N4xmmLUvT5PyTLu0jNaCqmQuAQwD+CHslYF2W22+/XQmpi8lTEx3dP1K5T3R0/0gnT02CrKuq\nOrp/pLgPc8vo/lF+uyaqo5GqSPRzkjiUqR2TiepgoBo5QKJlMJjdTjXaV3KdeBnlN6kUru3z3aft\nGi0aAM6qo8bm1uMWkTcA2AbwTUTW9pMATqrqD23bsB43qYvpuSne/+n348q1Kzc+W11exUPvfGju\nDaHIrDlLH1mCYv4ZEQiun7ZPSZNVVxuvM7dj/+e3cfG/zrdjNIrmunTZd5UByqxp1ZLtc4W1x2cJ\nOnWZiBwF8CUAb1HVJ0TkDIDLqvrPUuttANgAgLW1tdt3THeYkMDYJkQY7h/ie7/wvZnPisxTWXRu\nyyyRwynzPvH9EfDA/D5N05Y1kVVim1YNuDknpk8bQncEXSf0RArPAHhGVZ/Y+/tRAG9Kr6Sq26p6\nVFWP3nrrre6tJaQEtllsLr5wcS6oaBRLZKc5Fg3UZvmhrcc7ZP7cNG1ZnTPPZ7UjpkhmS5O++q6T\nK9yq+l0A3xaRV+99dAyR24SQVpMOKgrEuF5WmmPRQG3W3JHWsrIraxjM9hEYDCJLtg1sbWGufUl8\nM1s4v2ZxXLNKfg7AVESeAvAGAP+iuiYR4s5w/9D4uUDmRn8qdE68Xazn8ZExzp86j+unr+P8qfNO\n2TUmkYtF2GbFn3nHVqunLUtOq2bDx1rOukYkGyfhVtWv7rlBXqeqf0tV/2/VDSPEhTNvO4OVpZWZ\nz1aWVowBRSAS7zrSHLPmjsyy4ptwgfgQt88m3j7WMufXLA5neSeNMT03xebjm7hw6QLWDq1h69hW\nIRE17Wfz8c1CQUXiBjNCwsNZ3kllhBpFGHJgi8mV0fToz74PVKG13CwUbuJMSLGtugJhk6M/F6WW\nSJZbp+8dV9PQVUKcKZrTbKLowJY0odwtIVn0/GS6UYpBVwmphJDlXkNUIPR5A6izUNSi5yf7FsGi\nde4PhZs4E7LcawgftKu7pe5CUYuen+zTcS2KWyk0FG7iTMiAXwgftOsbQN2zyS96frJPx8UStcWg\ncBNnQgf8igxsSeL6BuDr4inrVln0jAufjmvR3UpFoXATLwqLbQWOTNc3AB8XTyi3StsH0lSJT8e1\n6G6lolC4SfVU5MjMewOILWdTnRKbi6dut0rTVBUYTHZcW1uR68N0jEV3KxWF6YCkehrIjzPV3hbI\njSHvtrTBUGmKXaCOtD2XY3Di44ig9biLQOEmM9gKOZsKTQeiaM65rb53H4fK19GfLnpOuw/M4ybt\nIpAj0ydoWCTnfHpuist/ennu89Xl1V5OlGyb62RnJ5zLhMHHaqBwk+oJ4Mj0DRoWyTnffHwTV69f\nnfv8pasvBYBWzfQeguVl+/9C5VIz+FgNFG5SPQHy43yDhkVyzm3W+MUXLtY6gKcurl2z/293F7j7\n7vIBSwYfq4HCTeqhZH6cr+ujSM65zRpfluVeZppkTYgQUzYBaNFz2quCwclFpyMh/ZAFrgBzcSoA\nxtnX06Id0/VME1PGhw0GE6uHwUniRocKRYQcbm/zlwOYs9KPv2wbcslsmhap0dI0ybztzU3g+HFg\naJ79bQYGE9sFLe5FpmO5WqFKuPpY7+vrwM6PTIG7NoDVRE74iwM88u566nuHIiunGoiE3JZp0tKv\nRK+gxU3c6Fiulmm4fZG6Iq7+8um5KXbetQ78zPuAq/uBHw4BFeD7I+inuyXaQHZBpzgEMZl0K5i4\nqCVh97msJCLnAfwAwDUAL7r2CqTlrK2ZTayO5GqlR0cmXR55QUiTxZ10fcT7xuE9pTtwEbgyAD71\nCHBu7BTYaxsu/XQc3uhA2GPuDSL29AHtbG9IfCzun1DVN1C0e4QpVwsAnn++E6ZL0boiLv5y076x\nugsc22y1BZqFrT9eWpq1WLtSIGuRS8LSVdI3fN4d41ytdHTq4sXWBimTFJ2RxyVV0LqPQxcy09nq\nnGnHF1s/fe1a62PTRjrm6QuKq3ArgC+IyJMislFlg0gJimSJjMfAwYPzn3fAdCk1I89TY+CB88BH\nrkc/n5pVYts+RofXMkU7na1y9yc2cMtPTAuJYboP/nu/ZO8UXPrrdE61aeRk+ra32Ye80KMyVTV3\nAfCKvZ9/FsDXALzVsM4GgLMAzq6trSlpgNFINZLs2WU0yt5OxLydiNfhJ09NdHT/SOU+0dH9I508\nNSl8Kq7HG2wNFPfhxjLYGuQedzJRHQxmT3UwiD4vs+/R/aOZ9W8sp0Zz+889t3Qbj0wUm+b2uJyP\nibzbXnS/ddH29vkC4Kw66LGq+qcDish9AJ5X1V+0rcN0wIYoWoUvQFqgqYzqYGVQaoYc1+P6pgi6\nnq7vvm0lYaECfOS6V0rdXBtPrQOHzSmMeOB8oduXdx26kC3akfFjTgQt6yoiBwAsqeoP9n7/PIB/\nrqqfs21D4W6Iok9agMLMoUc2VklVVWZt1wDfj8TVZ/9zbTy9BIjpWRXgPvNO846Xd9sbqMa70ITO\n4/5zAP6biHwNwO8C+C9Zok0apExFn/37b/4+HHoXlCgaKGyCqnyjpmwVXBkAj295739u3UvmjcXy\nucvx8uqILLQPueXkCreqfktVX7+3/EVV7WAiVM+wRYyKVPSJza6LickDXnjBu0mlAoU1U1XFujhb\nZbhvdGOgDj6zDZwbe+9/ro2PbwFXZxstLw6gXzDv1PV4Wal/Za9TmwObncfVGe6z3H777dV474k5\nIiOieuKE+/ajUbTNaKQ6HJojVMOhX7MKBgqbIn0ZQge0Quw/vY8TD84Gf3FkYrx1QLjzKXoepq/p\n6mr0tarqmncdeAQnKdxdw5Y5IpL/JJiepqzF5clKPNmTO4Y6+uiwtqwSZ0qqaNUiX7QNRZOI6sDW\ntr5kgFQBhbvP2HK4XJ5Yl6fJZ39dyMcq2cYmTjGdVnniwYmxDSdOtPfyZ31N29bJhCBE507h7jNZ\n4puXd+36NLnur80mX0zJNtZ9iiaXk/zTgdEtkrS82+Z+8LER2tZ2X0J17j7CzbKuXWM6Bd73vuj7\nkaZo4u7Skjm/a3kZePhhe3CzwnyxUCVcy7ax7pQ4a0rhtWVg6XqUXfL4FnBu3Oq0PJ9JGmI8M1Bb\nQ6h8d5Z17TPjMfCBD0TKkcQl3G9LE/jZn7UXscgaMl9RvpjvxMCZlGxj3Slx1vTJ5WtRHvfhnag2\n+JFpq9PybGVwsuhAlQUjTdRMoXBXRZW5UA8+CDzyiP9EfrZ0wQcfjH66FK9IUlFeXdGqf0ZKttG0\nuQhw552zn4UqLuWUPrm6C/nJTeMptK3IVTqzNG1vpOligagm8t0p3FVQx5RgRWtv2rYbj+3v3Ts7\n4fLGHQg6mMezjWnhw+umOH58VnBUIw9SfBlCviEYB/EY0EMX5k4h6JtKAExlV1WjW2CrZ97mtwgb\nTcxkT+GugrYXCra9DdieGhF7J1RB8ebgg3nGY0y3zmN97TqWLpzH+ubY2IdOz01xz6/fMyN89/z6\nPfjk09M5P3fydoZ8Q0iXnF0Ww1sQgJHhWgR9UwlAlguhCbGriiZmsqdwV0GbCwWb3gbuvhu45Zbo\n/d/kF8hSrQoIPTHwLVvruPt/Cnb+9j7ohwU771rHPffPl1o9+RsncfX61ZnPrl6/iot/6aRx3/Ht\nDD3cPzlF28Pvetj5WrSt7ECWC6EJsauSuiefoHBXQZuLPJjeBoBoyPvDD0fTfiefJlvWUYWdkMtE\nBy7EroOLL+4Agr0AH4DDO7j6NzZw8mOzyn3xhYvG/WBg/jy+nVUO9/e5Fm0rO9Anq7ptULiroM3f\n2CzB3d0FHnts1nTIc0ZWFIQNMTGwcfqxmNVdXHyD41uDZN/Oom8IrpfOdC1MhHxTCUGWVV1HGKjP\nULiroM3vgXlWf1rYszqhQE+fiyAXCbzluggOzf5/uN+eu7b/3nUM75gab2eRN4Qily5P6EO9qYTE\n5kJoexio7XAATtsJXSk+b2SEadSArQ01TsBQpN63dTDLHsN9I3xv8+a203NTvP/T78eVa1eM64ec\nGML30gUomd4qWOt7Hg7A6QtVvE9mjYywuXNsZlOAIKxrJoQ18Pb9HaufISu1blUGOPOO2XMdHxnj\noXc+FM0qY2D36i6O/+pmEK+Q76Xrm4Xa5jBQF6BwtxmXp7WIj3k8Bs6cmRfv5GQKLgR4+lwzIayB\nt0uwdmhJ1wGAG6l1o0MjPPQus+Uc+5MF5pEi1w5cCNKH+l66GUE/Mo2mMju9hJ13+Q+yaUOd7DaH\ngboAhbvN5JllRS1y0+QJQPS3jxoFePpcMyGM1rMCz68A0yOwmp+xEOtpxYsffhF6WjMDfHntSs5E\n42rxmoTS99LdEPQj02jI++GdG0PgfQbZtCUo2OYwUBegcLeZPLPMZpGfPJltUtlSAuPtXd+/Azx9\nTpkQ0ynGd21i+xO7GP4/wY35eAW4eADYuGtPvJMdWkmTMm8aspgLF7IPZxNKwO/S3RD6Y5vAavFB\nNravzPHj9Vvgdec+9wrXMoI+C8u6BiKvXqRrmdZ0jcm87fLKuYY+zVT96ZkJGFLXYHQKMyVP42V0\nKlHnNFCR6mS7lj84MpZWHQ6zD2crb7q87N+kyUQVp8V4/nKf2z1z+cq0pab3ogHW426IssWR5+aq\nOjE7tdhweHOfk0n09LsWPU4WkM4rllxlPW3fa5Rqq5w2C7ecRqXTwtj6A9vMb/HhsoSyiECO7h+Z\nO677zefnOlNdkcvV1lrgXYXC3QRlLT2XacXi/flOQZa2orO2j6dWqeKJLHKNUspntbg/OjSuH/It\nwiRUeYcL3Uf6zO1pm/dxZcXv62K7Fm2dfaerVCLcAJYBfAXAZ/PWXUjhLmvpuU4ZsrxsN5uy/pdu\nR9I6jS332Mo3iXrS2q/zGqW2mRyBDu6FXbhCW9w5ZmXe4fL62CL9SaZryaFtyVu+tFTscnVh8qOu\nUZVw/zyA/0DhtlDW0vOdVsy2lJ2IMOtpL2tSuVwjk7sodT6T21fskxKHNAUd9uVyuCyvVpVCl/eV\nWllR3bfP/Hne5arwxWZhCS7cAG4D8DiAv07htlCXxZ23lHV15D3tNsvd5VhFzNMi5xPK+ep4T10O\n14RroehXajgsvm9a3MWpQrgfBXA7gDtswg1gA8BZAGfX1tbqO9u2UIeP23Up8/S4PO2u5zyZzLpu\nDh6cd7C6pGDUoAZG90Ngs7LuYF7Rr5TL6dHHHZ6gwg3g7QAe3PvdKtzJZSEtblW/J9O0rslNYHvH\nzkoPKPO+euJE/lMdtzGrbbb27dsX/c90jRp6/7YG/O5wjBe0mKxEm7L9PrNKwhJauP8lgGcAnAfw\nXQC7ACZZ2yyscLviY65kuQ9sQufyrmsjlMumiDLUbHHHVrYpS+VGpkoPzMrJJMomSV/WpaX5zzt4\ner3BR7hzR06q6odU9TZVXQfwHgC/rap3+w/1ITfwqRhkG5342GPRs2biBz8oPvytjll6bMeosYBF\nskysjQsvPteLcdmbm8AVQ8HDl70MeOihzp/eQsIh703gWxrONDY4S2CvXDF3Ai5Dwesoz2Y7Rl7l\n/YCVkTInWYibeWgt2LjsJgs72b4qzz3XnmHnbSh81SlcTXOfha6SHEK4BPJcGmm/sKt7xhbRyvJb\n+4zgXF0tNqI0sMtC7jMPHc8b1FKEpgN5bc8Aafr6tAVw5GTLCfFNzUsZSBfD8Hl6bVGnrHa7+MaL\nDuKpQHkyfdsZg1piSmVBHpkoTo0Up7MH0IQiqy9ugzi2vWOpCwp3FwgRkk+n26UXl4JUvhkbtnZn\n5X+XCZaGbHvyNDyGjs9t69nvzjT/yERxb7HjlsH2VWmDZcvBPBEU7kXDZWhe1WaNbf8izQyVd8B1\n6HjZ5sysf2qUWSSqTH+e3Db2bCX3E+Iy2iolpEeKlqgjRoubwr1A5JktJjNxZcWeUx1vk/UEppUi\nPbhGJEpbLMtkMr9vl3HZFeFrIc5c+oyyrGU8aHmesxD1UrKOEV8T09cg7xzo446gcC8iLmZLWmiz\nknhdRkWaSs9ldQRFMSUiFwlyBqKIhXjDUs2wuMtYni4hBpd6KVl9dZkU/7xz4GAeCnc3CPVOnBx1\n6WO25KlE2f+XJXmOTVRoymlaYcs4w7dextfrM6dG0b66TB20RfNXF4HC3XZCvxMn63TndQZZw+iT\nT1jWU5pleoV4Ql2LbDSoBqX6XYtvvWqLO+nr9rGo8/pql2XR/NVFoHC3Hd8nNJT1eeyY+xNWZWm5\notfH45yLBh6bpGofd95+ioRJXJa8UAqJ6KZwL5KTy+edOJT1OZnk7yPvvdllOXCguuvjqERlUv2a\npsqskjx8wiTATRsifbuSoY68UAq5SfeEu49h5TLvpEkCWJ+qmj/ZoGk0RpHSckCJi5ZzzsvL+Uo0\nmejog8te8zL2gRBvGEUfw1Bf9UWne8Ldt7tbJCPD9oSEmJbbxdrOcoD6infV1y9nO+uEwo4zobcJ\nFws85BtG6BdfDq5xp3vC3be765ual/WEFLE+0/t2ndrbloDr4zax+bh9FaGIguxdK+uEwh2wuPNS\n4039l+/M73XSN5usSron3H27uyE7It8skqK+6SxLPG5H8lgnTpiLPJsG3JR1hSUzYZaX7YN69q57\n7oTCNePqxnC9denHwlYwK+QbRlFLvOytX6TQV/eEu28+7tAdkU/etqt17bpkdTamyRyMJmGJ62Gb\nkcck3onjTI5Elrecho4+uNyoaLu6MVw9UulbMvzoyCjcw4+OwpxDQ+LbN1nIo3vCrdqvrrWOb1yZ\npNoiFnfy3HxTE4vM7h5fK9sxlpfnr0kLn3QfN4brAJf05T3w5vnCVbh3oMM7wpx3Ey/Ek0nrxl1V\nTjeFu29U3RH5DmMbDmfbc+BA/jZp0Suampj35GcJbtZxylz3mgwFHzeGS19suiWAzpSKxalR9Lfl\nEnmfQ80hqLyvWVdDX3lQuBcB21M+dJgncTKZj3oB2RP5Zh0zS5DjbbJcKlnC7mNxu1KjZe5jcZua\nlVf+JW8Qa4hTqtvizvua0eKmcHeXMkPfbU+GSLYF6pOaaOsc4icvuf8sk87Hx+1KjUrkm6qXd+vS\n/8+7Hb4DWU3HN4Uyqqonppr9NaOPm8LdbSaT2UBk1gAan0kPsp4Qn9REW5A0bqdLumIspK5ZJa5k\nnX8FqhBq+L2pr3a5lWWCgaur5v53aclPUH08U1kdUltm7akCCnffcXnVN60TF43Ie9JtFqiPi8H1\nGDZ1iBWpCv9zljK02KQrGo92fZEoG++2HcfXM+Xi4w5R5r1tBBVuAC8B8LsAvgbgGwA+krcNhbti\nsvzbeRkfPoutdoqL6eR7rKWlm52KzR8eKqCYpwwlXCZVxjyLllV1DeaVKduadZwinqlkeMR2rJb2\nr4UJLdwC4ODe7ysAngDw5qxtKNwByFKAsk9YWRPKhSL55Fl56C5BV9/rm9WWgrusMuZZ1CJ29XNX\nZXGX6VDyqguHoC2ZyJW5SgAMAHwZwF/JWm9hhLuqO56nAHXkcOfNMOMSRTONrGxT55IVpC1wL8vG\nPF0uaREft+ssb64+btNntg5qMrG30eW65FndZSnT2YZ+/IMLN4BlAF8F8DyAf5W3fueF2+WOVGle\nueQ92zI2Qi1ZT3uW/zw9urOOTqYoZVUlRZl8Z9evUzpOe+xYsWHyWe0wDdJ1+cxE2b4x8C1ybl/e\nvqt4/Ku0uA8D+CKA1xr+twHgLICza2trxVvfNK53JEuQyna/LgrgMoAmuZgKIxd92l3EOGmx+4i3\nzSVic6GUdXba2lHAnHMVAZPouWxr+2qeOHFzfwFPJwhZbXLFtbJCyPblXa8qMkorzSoB8GEAH8xa\np9MWt+sdyXtHLfPNcrG4fQQ4fkp8LWDbt9dHhF2uVfqa2Uy8KkyvgE9g0WQf11GCLk2te7BM+vyL\ndEhF9x2Cou2rYjRp6ODkrQAO7/2+H8DvAHh71jadFm7XO+IigEWfFpMrJOm6KOJ+SO+/6Pt1loDa\njuvS3uXl/KfRt4NxvdYB33mLjn1yqcvhkn5edYA067xN/neTK6dNGZdFr1frLW4ArwPwFQBPAfg6\ngA/nbdNp4fZ5380TP19ByTJRkq6HInVK0sfJy/qwuSB8O42Q16oqc7LGtIK8UYFZAuKaft5ElkSW\nLzvpymkya8NGkevVKR+369Jp4fa5I3muBx9BcRG3eH8+qXZxfnSyjraLtZ1T89q7wzA5Kn2vVVPm\nZECy+p4iWSUh+68y1JG21zZan1Xiu3RauFX970gIQXGxZGMr2JZRMhzOmjemgGTZ/DFXizu9D1dz\nMY+2JN0WpOxXJSu80WTVvCxboq/V/EJD4W6CsoLiIqijkV0ATaMsyqTi2cykrI4ja5h6zfVBnKmp\nI0gepsgM7EmaDECasNUBa7pdXYPC3UXyRDY2y1yCpyHyp7PMpCzzymY+tk1tVGtzvYQ+TJs8Rnmx\n6o55shqFwt1FbGH5tAXrkiroEggsanGr5m9vyhCpojyr7Tq6WtA1dSZVHKYtHqM8+6Arot2G60nh\n7iou3548cyvvSVpdnR28c+DAvB88z0xyseZdUyJCiqSvKVpFMm5zh2mEPgQl2/IGQ+HuO1kO0ywx\nHQ7NRSZOnMiv7Z0+vm8eeB3q5ds5dNjibgtZaYBdsbbbcn8o3ItClnvF9C3MCmz6mhyuueAxdTwd\nvp1DR33cbcL2FexSvey2vBFRuBeFLHPHpBK+A3dcc6tdhv3VoV5FOocGsko6mMWYSdfPjRY3hbte\nXIQ46fbwzTTJMjnS7hpbnc+QeXB5FO0cWq48LW9e52nLGxGFe1HwDRL6uFayTA7TfpKzx8YCnWX9\nV0UTg6cqpOXN6w1t6Bwp3ItCkSBhMse7aAJuVofhUj62TVG5trwnW2h580hAfIR7CaS7jMfA9jYw\nGgEi9vUuXJjd5vz5aBtV8/qjUbTf8Th/f2kuXgSuXMlu984OsL4OLC1FP6fT7PWrxHYuWedYIy1v\nHmkICnfXiYX4+vVIcE2src1/ZnvyRaL92UTbtj8fRCLxVo1+bmw0J962cyl7joFoefOcmU7b01f3\nAQp3n9jaAgaD2c8Gg+jzNGUUwXQcH9KW/u4usLlZfH9l8LlmMTWqUJHmtY3pNOqb29JX9wJXn4rP\nQh93DiEjIel9uRY/zot65bUxOflhiKXKpFmXeqmu96OBaGEbAmdloJ/eDTA42WJCPvhlJ+OzKYKL\nqLsERX2KyvSJAAAK9klEQVSWIk9xiBIBvlCFvGnLAJe2Q+FuM6Ee/KwBNWUHnNjauLwcbRfS0i4q\npK6CHFpoqULesK9zg8LdZkI9+EUq17uKne8IS99lOCz/7u+qBqGFlirkDXPR3fARbgYn6yZEmsB0\nGqXduRwjGUg7fjwKBCZJBwan02jdMuQFLs+cuZkJk5fBYsM1Ty50WkYfooU1k85azcs2JQ64KrzP\nQos7gxDmR9YAmGRZNh9ftG3ouu8S1+K2tfHgwTCRNlfL1/N6OwUCHVbqekCR1A8Cz/L+SgBfBPBN\nAN8AcDJvm94Ld9mnsuz2Wa6MZFk239okZZe8TmN11V7TxBebIJuyahyvd6hXepf9UNhJmtDC/aMA\n3rT3+0sB/AGA12Rt02vhboPDLqs8a5KqfdWmJUlanWx++aL+YVMqZIl7E8p9nbefNnyFSPsIKtxz\nGwCfBvBTWev0WribDk7Z6mD7ZFTE2SGhRTvvGlSdkVHy3oRqXt5+mv4KkXbiI9xeUSgRWQfwRgBP\nlPWtd5Ymi0fEQ9DSgcnh0BztsQXSHn44e4i8C8vLs3/Hw9jjkYSm0YVVj98ueW9CNS9vP6w/Qkrj\nqvAADgJ4EsDPWP6/AeAsgLNra2u19VK106S5FHqigKzyrLFlbrOu45S+2JRM/m9lxTyPZUlXRq5j\nuOS9qcvHTYubmEBoVwmAFQC/CeDnXdbvtaukSQdl0Xf5PPG2jZ7My15R9QuAxvsvEpVzjfiVvDeh\ngoa+/SV93CSocAMQAL8K4AHXnfZauFWbSwkoanEXmU8yL40wPqaPr7yML9sn/a8D6RqTierwjoni\n1EhxWnT40ZFOnmpnW0k9+Ai3ROvbEZG/BuB3AJwDcH3v43tV9THbNkePHtWzZ88W890QO7GPOzmI\nZnUVeOlLgeeei5yoW1uzvu719cj3nGY0iga/mLBtEzMY3PSp563resw8lpYiqU4jEvnrO8b03BQb\nn9nA7tWb93KwMsD2XdsYH+HIlEVERJ5U1aNO6+YJdxEo3BUynUYjHS9cAF7+cuDyZeDq1Zv/T4oq\nUEzwbNsAUVDy2rVIhOPRgunOxIQI8MgjxYfLFemAWsz6A+vYuTR/PqNDI5w/db7+BpHG8RFuDnnv\nGsmJEw4enBVtYH4Iuy3F4eUvtx/Dto1IJNrAzaLKwM3xzFmolhvj3LOh5hcumVNIbJ8TkoTC3WVc\n8sq2tiJ3SprLl+2V7E0iKTJvhe/uAidP3uxMJhP7FGplUg+B3hW8WDtk7hxtnxOShMLdZVwSj8fj\nyAee5upV+6wzJpG0uU4uXrzZAYzHwAc+MC/eoSzj5NtG0eJULWHr2BYGK7Od42BlgK1j3XyDIPVC\n4e4yru6D554zb5814mM8jvazthatlx5wkyTZATz4YOTL7ollXBXjI2Ns37WN0aERBILRoREDk8Qd\n1/QTn6X36YBtIi8PO2viA980wqylioRnQhYIcCKFBcdFdPNyuX0rCxYpf8eRKITcwEe4mQ7YR2yp\nc8vLkX/YlO+dJisl0EacmmfKN0+nKWa1s6MpfoSUgXnci06IwSo+A2vS+3cV5J4NqiGkDMzjXnSK\nlLlLV/O78878Kchs+29qWjFCFgQKdx/xHawSuzZ2diILeGcnKv16/Lh7/vXq6s39uwpyzwbVEFIX\nFO4+4jtYZXPTPInwJz9582/bwJqYpMvDVZB7NqiGkLqgj7sLJOuTuAQWfXENRMajJ+N6JWmSPuyq\n20xIz2Bwsk+4ZmiUwbfC34ULDCoSEhgGJ/uEzY1hG65eBJNrw0ZsQZvoS1DRNO0aIS2Cwt12Qk1Q\nmCVGJl/zcGjeT+z26GtQ0RSo3digeJNWQeFuOyGs2yJi9O5328W5z0HFOt5wCCkJfdxtJ4SPO29A\njO0Yx48Djz22WAFGDgoiDeHj495XdWNISWKhLJOhkedusVmZjz22eEPP19bMnVxf/PekF9BV0gXK\n1qHOc7eE8qP70sYgYJ/996Q3ULj7znQKPP/8/OdJMWoiSyTL796koPfZf0/6Q175QAAPAfgTAF93\nLTnIsq4twVbedTicL8Eaqryqa31tW9nY4ZClXslCgpD1uAG8FcCbKNwdxCaOpgkUQkxo4NMBiPjV\n+86a9IGQHuAj3E5ZJSKyDuCzqvpaFyueWSUtoe4MCZ/62r5lY5nVQXoOR06SiLp91z5BTlsQMGvg\nDyEEQEDhFpENETkrImefffbZULslZag7Q8Kno7AFAc+cYVYHIXm4+FMArIM+7m5S52S8oYKcnECY\nLCCgj5s0Bsu5ElKIoD5uEfk4gP8O4NUi8oyI/N2yDSQ9puxgoSK0cSAPIRWSK9yq+l5V/VFVXVHV\n21T1l+to2EJSpwD1RexYzY8sICwy1RbqmDChiWNVjU8KIiEthjPgdJE6BahPYsdqfqQnMI+7i9RZ\n6KmpolJV0PfZeAgxQOFuC3UKUJ/EjtX8yAJC4W4LdQpQn8SO1fzIAkLhbgt1ClDfxK6JFERCGoTB\nSUIIaQEMThJCSI+hcBNCSMegcBNCSMegcBNCSMegcBNCSMeoJKtERJ4F8EMA3wu+83ZxC3iOfWER\nzpPn2G5Gqnqry4qVCDcAiMhZ19SWrsJz7A+LcJ48x/5AVwkhhHQMCjchhHSMKoV7u8J9twWeY39Y\nhPPkOfaEynzchBBCqoGuEkII6RiVCreI/BsR+T0ReUpEfk1EDld5vDoRkZ8Wkd8XkT8UkX/SdHtC\nIyKvFJEvisg3ReQbInKy6TZVhYgsi8hXROSzTbelKkTksIg8uvc8Pi0iP950m0IjIv9o77v6dRH5\nuIi8pOk2VUXVFvfnAbxWVV8H4A8AfKji49WCiCwD+PcA3gbgNQDeKyKvabZVwXkRwD9W1dcAeDOA\nv9/Dc4w5CeDpphtRMWcAfE5VfwzA69Gz8xWRVwD4hwCOquprASwDeE+zraqOSoVbVX9LVV/c+/NL\nAG6r8ng18pcB/KGqfktVrwD4BIB3NtymoKjqH6vql/d+/wGiB/0VzbYqPCJyG4C/CeBjTbelKkTk\nEIC3AvhlAFDVK6r6/WZbVQn7AOwXkX0ABgD+d8PtqYw6fdzvB/AbNR6vSl4B4NuJv59BD0UtRkTW\nAbwRwBPNtqQSHgDwCwD6PLPwqwA8C+BX9lxCHxORA003KiSq+h0AvwjgAoA/BnBJVX+r2VZVR2nh\nFpEv7PmU0ss7E+tsInr1npY9HqkXETkI4D8BOKWql5tuT0hE5O0A/kRVn2y6LRWzD8CbAPySqr4R\nUTmKXsVlRORliN56XwXgzwM4ICJ3N9uq6thXdgeq+pNZ/xeRvwPg7QCOaX9yD78D4JWJv2/b+6xX\niMgKItGequqnmm5PBbwFwDtE5E4ALwHwIyIyUdW+PfDPAHhGVeM3pkfRM+EG8JMA/khVnwUAEfkU\ngL8KYNJoqyqi6qySn0b0GvoOVd2t8lg18z8A/AUReZWIrCIKgvznhtsUFBERRD7Rp1X13zbdnipQ\n1Q+p6m2quo7oHv52D0UbqvpdAN8WkVfvfXQMwDcbbFIVXADwZhEZ7H13j6FnAdgkpS3uHP4dgD8D\n4PPRtcSXVPUDFR+zclT1RRH5BwB+E1H0+iFV/UbDzQrNWwC8D8A5Efnq3mf3qupjDbaJFOfnAEz3\nDI1vAbin4fYERVWfEJFHAXwZkVv2K+jxKEqOnCSEkI7BkZOEENIxKNyEENIxKNyEENIxKNyEENIx\nKNyEENIxKNyEENIxKNyEENIxKNyEENIx/j/qQxVBlPXeoQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f4db876f400>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADGtJREFUeJzt3V+InfWdx/HPp0mKNts2Fzk7ExzT2QvphdJVGVJbS2mz\nWPwT9MaLXFjBm5AiiwtbZN2LgvfL4tqCQ7AUxRaRthYJiZCiZetFlBmNsTa5CCXFBCcZC8amFhf1\nsxfzFKaTSc5zZk58zvn6fsFhnvOc35zz5UHfefL4zOgkAgDU8pmuBwAADB9xB4CCiDsAFETcAaAg\n4g4ABRF3ACiIuANAQcQdAAoi7gBQ0MauPnjr1q2Znp7u6uMBYCzNz8+/k6TXb11ncZ+entbc3FxX\nHw8AY8n2H9us47IMABRE3AGgIOIOAAURdwAoiLgDQEHEHeNvclKyL3xMTnY9GdCZVnG3fdL2G7aP\n2L7g/kUvedT2CdtHbd84/FGBizhzZrD9wKfAIPe5fzvJOxd57TZJ1zSPr0p6rPkKAOjAsC7L3CXp\nySw5LGmL7W1Dem8AwIDaxj2Sfm173vaeVV6/StJby56favYBADrQ9rLMN5Kctv2Pkg7ZPp7kfwf9\nsOYPhj2StH379kG/HQDQUqsz9ySnm69nJT0raceKJaclXb3s+VSzb+X77Esyk2Sm1+v7e2+AdiYm\nBtsPfAr0jbvtzbY//7dtSd+R9LsVy56TdG9z18xNks4leXvo0wKrWViQkgsfCwtdTwZ0ps1lmQlJ\nz9r+2/qfJXne9l5JSjIr6YCk2yWdkPS+pPsuz7gAgDb6xj3JHyT98yr7Z5dtR9L9wx0NALBW/IQq\nABRE3AGgIOIOAAURdwAoiLgDQEHEHQAKIu4AUBBxB4CCiDsAFETcAaAg4g4ABRF3ACiIuANAQcQd\nAAoi7gBQEHEHgIKIOwAURNwBoCDiDgAFEXcAKIi4A0BBxB0ACiLuAFAQcQeAgog7ABRE3AGgIOIO\nAAURdwAoiLgDQEHEHQAKIu4AUBBxB4CCiDsAFETcAaAg4g4ABRF3ACiIuANAQcQdAAoi7gBQUOu4\n295g+zXb+1d57Vu2z9k+0jx+MNwxAQCD2DjA2gckHZP0hYu8/tsku9Y/EgBgvVqdudueknSHpMcv\n7zgAgGFoe1nmEUkPSvr4Emu+bvuo7YO2r11tge09tudszy0uLg46KwCgpb5xt71L0tkk85dY9qqk\n7Um+IumHkn612qIk+5LMJJnp9XprGhgA0F+bM/ebJd1p+6SkpyXttP3U8gVJ3ktyvtk+IGmT7a3D\nHhYA0E7fuCd5KMlUkmlJuyW9kOSe5WtsT9p2s72jed8/XYZ5AQAtDHK3zN+xvVeSksxKulvS92x/\nKOmvknYnyXBGBAAMyl01eGZmJnNzc518NgCMK9vzSWb6reMnVAGgIOIOAAURdwAoiLgDQEHEHQAK\nIu4AUBBxB4CCiDsAFETcAaAg4g4ABRF3ACiIuANAQcQdAAoi7gBQEHEHgIKIOwAURNwBoCDiDgAF\nEXcAKIi4A0BBxB0ACiLuAFAQcQeAgog7ABRE3AGgIOIOAAURdwAoiLgDQEHEHQAKIu4AUBBxB4CC\niDsAFETcAaAg4g4ABRF3ACiIuANAQcQdAAoi7gBQEHEHgIJax932Btuv2d6/ymu2/ajtE7aP2r5x\nuGMCAAYxyJn7A5KOXeS12yRd0zz2SHpsnXMBANahVdxtT0m6Q9LjF1lyl6Qns+SwpC22tw1pRgDA\ngNqeuT8i6UFJH1/k9askvbXs+alm39+xvcf2nO25xcXFgQYFALTXN+62d0k6m2R+vR+WZF+SmSQz\nvV5vvW8HALiINmfuN0u60/ZJSU9L2mn7qRVrTku6etnzqWYfAKADfeOe5KEkU0mmJe2W9EKSe1Ys\ne07Svc1dMzdJOpfk7eGPCwBoY+Nav9H2XklKMivpgKTbJZ2Q9L6k+4YyHQBgTQaKe5LfSPpNsz27\nbH8k3T/MwQAAa8dPqAJAQcQdAAoi7gBQEHEHgIKIOwAUtOZbIQEA7U3+16TO/OXMBfsnNk9o4fsL\nQ/88ztwB4BOwWtgvtX+9iDsAFETcAaAg4g4ABRF3ACiIuAPAJ2Bi88RA+9eLWyEB4BNwOW53vBTO\n3AGgIOIOAAURdwAoiLgDQEHEHQAKIu4AUBBxB4CCiDsAFETcAaAg4g4ABRF3ACiIuANAQcQdAAoi\n7gBQEHEHgIKIOwAURNwBoCDiDgAFEXcAKIi4A0BBxB0ACiLuAFAQcQeAgog7ABRE3AGgoL5xt32F\n7Vdsv277TdsPr7LmW7bP2T7SPH5wecYFALSxscWaDyTtTHLe9iZJL9k+mOTwinW/TbJr+CMCAAbV\nN+5JIul883RT88jlHAoAsD6trrnb3mD7iKSzkg4leXmVZV+3fdT2QdvXDnVKAMBAWsU9yUdJrpc0\nJWmH7etWLHlV0vYkX5H0Q0m/Wu19bO+xPWd7bnFxcT1zAwAuYaC7ZZK8K+lFSbeu2P9ekvPN9gFJ\nm2xvXeX79yWZSTLT6/XWMTYA4FLa3C3Ts72l2b5S0i2Sjq9YM2nbzfaO5n3/NPxxAQBttLlbZpuk\nJ2xv0FK0n0my3/ZeSUoyK+luSd+z/aGkv0ra3fyHWABAB9rcLXNU0g2r7J9dtv0jST8a7mgAgLXi\nJ1QBoCDiDgAFEXcAKIi4A0BBxB0ACiLuAFAQcQeAgog7ABRE3AGgIOIOAAURdwAoiLgDQEHEHQAK\nIu4AUBBxB4CCiDsAFETcAaAg4g4ABRF3ACiIuANAQcQdAAoi7gBQEHEHgIKIOwAURNwBoCDiDgAF\nEXcAKIi4A0BBxB0ACiLuAFAQcQeAgog7ABRE3AGgIOIOAAURdwAoiLgDQEHEHQAKIu4AUBBxB4CC\n+sbd9hW2X7H9uu03bT+8yhrbftT2CdtHbd94ecYFALSxscWaDyTtTHLe9iZJL9k+mOTwsjW3Sbqm\neXxV0mPNVwBAB/qeuWfJ+ebppuaRFcvukvRks/awpC22tw13VABAW62uudveYPuIpLOSDiV5ecWS\nqyS9tez5qWYfAKADreKe5KMk10uakrTD9nVr+TDbe2zP2Z5bXFxcy1sAAFoY6G6ZJO9KelHSrSte\nOi3p6mXPp5p9K79/X5KZJDO9Xm/QWQEALbW5W6Zne0uzfaWkWyQdX7HsOUn3NnfN3CTpXJK3hzno\n5KRkX/iYnBzmpwBADW3ultkm6QnbG7T0h8EzSfbb3itJSWYlHZB0u6QTkt6XdN+wBz1zZrD9APBp\n1jfuSY5KumGV/bPLtiPp/uGOBgBYK35CFQAKIu4AUBBxB4CCxibuExOD7QeAT7M2d8uMhIWFricA\ngPExNmfuAID2iDsAFETcAaAg4g4ABRF3ACjIS785oIMPthcl/XEN37pV0jtDHqcijlN7HKt2OE7t\nXO7j9KUkfX+tbmdxXyvbc0lmup5j1HGc2uNYtcNxamdUjhOXZQCgIOIOAAWNY9z3dT3AmOA4tcex\naofj1M5IHKexu+YOAOhvHM/cAQB9jE3cbV9t+0Xbv7f9pu0Hup5pFNm+wvYrtl9vjtPDXc80ymxv\nsP2a7f1dzzKqbJ+0/YbtI7bnup5nlNneYvvnto/bPmb7a13NMja/FVLSh5L+Pcmrtj8vad72oSS/\n73qwEfOBpJ1JztveJOkl2weTHO56sBH1gKRjkr7Q9SAj7ttJuMe9v/+R9HySu21/VtLnuhpkbM7c\nk7yd5NVm+89a+hfyqm6nGj1Zcr55uql58B9WVmF7StIdkh7vehaMP9tflPRNST+WpCT/l+TdruYZ\nm7gvZ3taS//T7pe7nWQ0NZcajkg6K+lQEo7T6h6R9KCkj7seZMRF0q9tz9ve0/UwI+yfJC1K+klz\nqe9x25u7Gmbs4m77HyT9QtK/JXmv63lGUZKPklwvaUrSDtvXdT3TqLG9S9LZJPNdzzIGvtH883Sb\npPttf7PrgUbURkk3SnosyQ2S/iLpP7oaZqzi3lxD/oWknyb5ZdfzjLrmr4QvSrq161lG0M2S7rR9\nUtLTknbafqrbkUZTktPN17OSnpW0o9uJRtYpSaeW/U3551qKfSfGJu62raVrWceS/HfX84wq2z3b\nW5rtKyXdIul4t1ONniQPJZlKMi1pt6QXktzT8Vgjx/bm5gYGNZcYviPpd91ONZqSLEh6y/aXm13/\nIqmzGz7G6W6ZmyV9V9IbzfVkSfrPJAc6nGkUbZP0hO0NWvrD+5kk3OaHtZqQ9OzSuZU2SvpZkue7\nHWmk/auknzZ3yvxB0n1dDcJPqAJAQWNzWQYA0B5xB4CCiDsAFETcAaAg4g4ABRF3ACiIuANAQcQd\nAAr6f1g3CT2eCeGJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f4db875ce48>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "finite_mixture_model()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
